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We characterize empirically the radial diffusion of stars in the plane of a typical barred disk galaxy by calculating the local spatial 
diffusion coefficient and diffusion time- scale for bulge-disk-halo A/^-body self-consistent systems which initially differ in the Safronov- 
Toomre-2r parameter. We find different diffusion scenarios that depend on the bar strength and on the degree of instability of the 
disk. Marginally stable disks, with 2r ~ 1, have two families of bar orbits with different values of angular momentum and energy, 
which determine a large diffusion in the corotation region. In hot disks, Qt > ^, stellar diffusion is reduced with respect to the case 
of marginally stable disks. In cold models, we find that spatial diffusion is not constant in time and strongly depends on the activity 
of the bar, which can move stars all over the disk recurrently. We conclude that to realistically study the impact of radial migration on 
the chemical evolution modeling of the Milky Way the role of the bar has to be taken into account. 

Key words. Galaxies: kinematics and dynamics. Galaxies: stellar content. Galaxies: spiral. Galaxy: disk, bulge. Methods: numerical 
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1. Introduction 

Disk galaxies are highly nonlinear systems which are driven by 
external forcing (satellites, interaction with nearby galaxies), in- 
ternal instabilities (related to the formation of internal structures 
such as spiral arms, rings, bars) or a combination of both. Chaos 
and complexity, which are two different aspects of nonlinear re- 
sponse, dominate the dynamics of galactic systems. Disk galax- 
ies are 'complex' in the sense that they are made of many com- 
ponents (stars, gas, dark matter in the bar-bulge, disk and halo 
components) whose interactions can give rise to spontaneous 
self-organization and to the emergence of coherent, collective 
phenomena. Examples of emergent behavior in disk galaxies are 
the formation of a central bar or the onset of spiral arms and 
warps, which can develop in dynamically cold disks (Revaz & 
Pfenniger 2004). Disk galaxies are not only complex systems, 
but even display chaotic behavior, since tiny differences in initial 
orbits of stars can exponentially blow up, as indicated by numer- 
ical simulations. Since the chaotic orbits are more sensitive to 
perturbations than regular (periodic) ones, external/internal forc- 
ing are more effective on this chaotic component. In this paper, 
we investigate chaotic and complex phenomena related to the 
formation of a central bar which give rise to the diffusion of the 
stellar component in the disk. 

The bar has a time-dependent activity, with a pattern speed 
which typically decreases in isolated galaxies (Sellwood 1981). 
However, the system can be cooled or heated by energy dissipa- 
tion or infall of gas, or by forming stars on low- velocity dis- 
persion orbits, with the net effect of impacting the amplitude 
of spiral waves and the strength of the bar, or even destroy- 
ing it. In this way bars (and spiral waves) can be seen as re- 
current patterns which can be rebuilt during their long history 
until the present configuration at redshift z = (Bournaud & 
Combes"2002 ). Under the action of these non-axisymmetric pat- 
terns, stars move in the disk and gradually increase their velocity 
dispersion, as suggested by observations in the Solar neighbor- 
hood (see Holmberg et al. <2009i and references therein) and in 



external galaxies (Gerssen et al. 2000', Shapiro et al. 12003b . The 
origin and the amount of disk heating are still open to debate. 

First attempts to explain such a heating process in disk galax- 
ies were made by empirically modeling the observed increase 
of the stellar velocity dispersion with age in the solar neighbor- 
hood. Wielen (11977b suggested a diffusion mechanism in veloc- 
ity space, which gives rise to typical relaxation times for young 
disk stars of the order of the period of revolution and to a devi- 
ation of stellar positions of 1.5 kpc in 200 Myr. The result was 
obtained without making detailed assumptions on the underly- 
ing local acceleration process responsible for the diffusion of 
stellar orbits. Global acceleration processes, such as the gravi- 
tational field of stationary density waves or of central bars with 
constant pattern speed, were ruled out since their contribution to 
the velocity dispersion of old stars was found to be negligible 
and concentrated in particular resonance regions (Wielen [19771 
Binney & Tremaine 2008, p. 693). In isolated galaxies, differ- 
ent local accelerating mechanisms have been investigated, such 
as the gravitational encounters between stars and giant molecu- 
lar clouds (Spitzer & Schwarschild 1951. ,1953^ Lacey .1984.) , 
secular heating produced by transient spiral arms (Barbanis 
& Woltjer 1967; Carlberg & Sellwood T985; Fuchs ^20®or 
the combination of the two processes (Binney & Lacey 119881 
Jenkins & Binnev 19901) . Another heating mechanism was sug- 
gested by Minchev and Quillen ('2006), who showed that the 
stellar velocity dispersion can increase with time due to the non- 
linear coupling between two spiral density waves. 

Such local acceleration mechanisms suggest the existence of 
a significant component of the galactic gravitational field with 
a rather chaotic behavior. Pfenniger (T986) investigated the re- 
lation between diffusion and chaotic orbits. The latter typically 
react promptly to small perturbations. He pointed out that the ef- 
fect of perturbations on regular orbits, such as epicyclic orbits, 
underestimates strongly the stellar diffusion rate when a stellar 
system becomes non-integrable, as for example in the presence 
of a central bar. As the central bar develops, reaches its maximal 
amplitude and then settles down to an almost steady state, its 
gravitational potential changes in time. In time-dependent po- 
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tentials, the number of chaotic orbits typically decreases while 
the system secularly evolves toward a quasi-steady state through 
collective effects. The system is then ready again to respond 
(mainly through the remaining irregular orbits) to external per- 
turbations, such as new infall of gas, and to recurrently restore 
a strong bar. The bar is thus able to perturb orbits of stars born 
or passing through its region, which can visit at later times the 
Solar neighborhood (Raboud et al. 119981) . Indeed, most of the 
observational signatures of radial mixing reported in the litera- 
ture (Grenon rT972l 119991 Castro et slI^TWT) point to stars com- 
ing from a region next to the bulge/bar intersection, suggesting 
the bar to be a key player in the radial migration process. 

The subject of radial migration of stars was revived when 
a large scatter in the observed age-metallicity relation (AMR) 
was reported by Edvardsson et al. (119931) , later confirmed by 
the la rger Geneva-Cop enhage n Survey sample (Nordstrom et al. 
[20041 Holmberg et al. I^QTl (2009'; see also Casagrande et al. 
12011b . It must be said that even though the AMR has been ex- 
tensively studied in the solar vicinity, the results are still con- 
troversial due to the large uncertainties in stellar ages (see Pont 
& Eyer 2004). For instance, using a sample of stars for which 
it was possible to obtain chromospheric ages, Rocha-Pinto et al. 
(12000) have reported a much tighter AMR. 

The mechanisms driving radial diffusion and heating are still 
hotly debated, and in many cases the role of the bar is not taken 
into account. Assuming the whole scatter seen in Edvardsson 
et al. ([T993I) data was real, Sellwood & Binney (l2QQ2b po inted 
out that the radial excursion predicted by Wielen (11977) was 
not sufficient to explain the weakness of the AMR in the so- 
lar neighborhood. In order to explain both the large scatter in 
the AMR and the evidence that even old disk stars today have 
nearly circular orbits, Sellwood & Binney (12002b suggested a 
new mechanism based on the resonant scattering of stars under 
the effect of transient spiral waves. In this process, a star ini- 
tially on a nearly circular orbit resonates with a rotating wave 
and changes its angular momentum. If the duration of the peak 
amplitude of the perturbing potential is less than the period of 
the 'horseshoe' orbits, i.e. orbits of particles trapped at the coro- 
tation radius of the spiral wave, the star can escape from the 
potential well without changing its eccentricity. The net effect of 
this scattering mechanism is that stars migrate radially without 
heating the disk. In other words, the overall distribution of an- 
gular momentum is preserved, except near the corotation region 
of the transient spiral wave, where stars can have large changes 
of their angular momenta. Haywood (12008b estimated upper val- 
ues for the migration rate from 1.5 to 3.7 kpc/Gyr, which agree 
with the values in Lepine et al. (12003b for the radial wandering 
due to the scattering mechanism assumed by Sellwood & Binney 
(I2QQ21) . 

Radial diffusion of stars (and gas) could have important im- 
plications for the interpretation of key observational constraints 
for the formation of the Galaxy, such as the AMR, metallicity 
distributions, or the metallicity gradients, since old, probably 
more metal rich stars that formed at small galactocentric radii, as 
well as young metal-poor stars formed at large radii are enabled 
to appear in Solar-neighborhood samples (e.g. Haywood 2008 ). 
Due to the lack of detailed information on the processes driving 
stellar radial migration, models of the Galactic chemical evolu- 
tion have evaluated past history of the solar neighborhood and 
the formation and evolution of the abundance gradients assum- 
ing that the Galaxy can be divided into concentric wide (~1- 
2 kpc) cylindrical annuli, which evolve independently (van den 
BerghflMl S chmid t gggH Pagel fTWTl Chiappini et al. 1 19971 
Chiappini et al. 1200 lb . Schonrich & Binney ([2009aJ explored the 



consequences of mass exchanges between annuli by taking into 
account the effect of the resonant scattering of stars described 
before. This approach appears to be successful to replicate many 
properties of the thick disk in the Solar neighborhood without re- 
quiring any merger or tidal event (Schonrich & Binnev l2009bb . 
High resolution cosmological simulations (Roskar et al. I2008[ 
Loebman et al. 2010) give support to the view that such scatter- 
ing mechanism determines a significant migration in the stellar 
disk. However, the strong mixing driven by bar resonances was 
not taken into account (see below), casting thus doubts on some 
of the conclusions in the papers quoted above. 

The Milky Way (MW) is a barred galaxy and it is clear 
that the process above, not accounting for the existence of the 
bar, is probably just one of the processes at play among others. 
Indeed, the role of resonant couplings between bars and spirals 
(Tagger et al. |1987 ) in the distribution of energy and angular mo- 
mentum in disk galaxies could also play a major role. Recently, 
Minchev and collaborators (Minchev & Famaey 2 10^ Minchev 
et al. 12011b have further analyzed this mixing mechanism find- 
ing that resonances between the bar and the spiral arms can 
act much more efficiently than transient spiral structures, dra- 
matically reducing the predicted mixing time- scales. Moreover, 
while for the Sellwood & Binney (2002) mechanism to work 
short-lived transient spirals are required, in barred galaxies, such 
as the MW, spirals are most likely coupled with the bar as shown 
by Sparke & Sellwood 1987, and thus longer lived (Binney & 
Tremaine 2008, Quillen et al. 12010b . As a consequence the ra- 
dial migration process in the MW could have been different than 
currently predicted. 

In order to include the effect of radial migration in chemical 
evolution models and to gain a global (chemical and kinematic) 
understanding of the processes at play in the galactic disks, many 
dynamical aspects need to be further investigated and in partic- 
ular the role of the bar, that is the strongest non-axisymmetric 
component in disk galaxies. In this paper, we present A/^-body 
simulations of barred spiral galaxies, and study how disks with 
different degrees of stability, ranging from marginally stable 
disks with Safronov-Toomre parameter ~ 1 to hot disks with 
Qt > I, respond to the presence of bar patterns. Our aim is to 
estimate the time and length scales of stellar diffusion in the ra- 
dial direction and to relate these quantities to the strength of the 
bar and to the number of hot particles in the disk, i.e. generally 
chaotic particles which are susceptible to cross the corotation 
barrier and to explore all space, being characterized by values 
of the Jacobi integral H larger than the value at the Lagrangian 
points, H > H{Li2) (Sparke & Sellwood 1987; Pfenniger & 
Friedli 11991b . We investigate how these characteristic scales 
evolve in time, and depend on the activity of the bar. We consider 
different scenarios of diffusion, and discuss their implications for 
chemical evolution constraints in our Galaxy. 

The paper is organized as follows. In Sect. 2 we describe 
the simulations and the relevant parameters. In Sect. 3 we solve 
the diffusion equation in axisymmetric systems, we define the 
diffusion coefficient, the diffusion time- scale and the diffusion 
length-scale, and the methods used to estimate these quantities 
from the simulation results. In Sect. 4 we present our results. 
In Sect. 5 we discuss the implications for chemical evolution 
models of the MW and we summarise our findings. 

2. A^-body simulations 

We have run self-consistent A/^-body simulations starting from a 
bar-unstable axisymmetric model. We have analyzed initial con- 
figurations with disk, bulge and dark halo components which 
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differ on the initial value of the Safro nov-Too mre pa rameter 
Qj = crrK/(336Gi:) (Safronov I960'; Toomre T96T), where 
CTf is the radial velocity dispersion of the disk component, G 
is the gravitational constant, S is the disk surface density, k is the 
epicycle frequency defined by = RdD?^/dR + 4Q^, where Q is 
the circular frequency related to the global gravitational potential 
0(7?, z, t) in the disk plane z = Ohy Q.^ = (l/R) d^/dR. 

The initial mass distribution in our simulations corresponds 
to a superposition of a pair of axisymmetric Miyamoto-Nagai 
disks of mass Mb, Md, horizontal scales Ab -\- B, Ad -\- B, and 
identical scale-height B, 



-GMi 



(1) 



The first component represents the bulge (B), while the second 
the disk (D) (Pfenniger & Friedli 1991). The parameters have 
been set to Afi = 0.07 kpc. A/) = 1.5kpc, 5 = 0.5 kpc, Mg/M/) = 
3/17. The initial particle positions and velocities are found by a 
Monte-Carlo draw following the density law corresponding to 
Eq. ([T]), truncated to a spheroid of semi-axes R = 30 kpc, z = 
10 kpc. The number of particles in the disk-bulge component is 
= 4 • 10^ and the total mass is Mbd = 4.2 • 10^^ M©. 
In order to progressively heat the disk, we have added to this 
bulge-disk component an oblate pseudo-isothermal halo with the 
following density distribution (except in models ml and m2): 



Ph(R,z) = 



Ph 



1+ 7^2/^2 +^2/^2 



(2) 



The number of particles in this halo component is Nh = 2 • 
10^, which is a value in the range suggested in Dubinski et al. 
2009 ) in order to obtain convergent behavior in studies of bar 
formation and evolution. The length-scales are Rh = 7.5 kpc 
and Zh = 3.5 kpc. The density distribution has been truncated to 
R = 30 kpc, z = 15 kpc. We set the total mass in the dark halo 
Mh to be four times the total mass in the bulge-disk component 
Mbd, except in the model m3, where Mh = 2Mbd (see Table [T]). 
The effect of adding the halo component is that the bar becomes 
progressively smaller and with higher pattern speed, the disk is 
hotter and less sensitive to the bar perturbations. 

We then impose the equilibrium of the first and second mo- 
ments of velocities by solving Jeans' equations (see e.g., Binney 
& Tremaine 2008 ) with a constant Qt. The resulting distribu- 
tion is relaxed for a couple of rotations until ripples spread- 
ing through the disk from the center disappear. We use this 
as the initial condition for the A/^-body simulations performed 
by using the Gadget- 2 free source code (Springel et al. 120011 
Springel 2005 ). 

The initial Gadget- 2 configurations considered in this work 
differ on the values of the Safronov-Toomre parameter, which 
ranges from 2^ ~ 5 at two scale lengths from the center for hot 
disks to 2r ~ 1 for marginally stable disks. These Qt values 
are listed in Table [H along with the initial values of the radial 
and vertical velocity dispersions at two scale lengths from the 
center. We have considered these initial values in order to inves- 
tigate how the radial diffusion depends on the disk sensitivity to 
perturbations. Thus, we have considered two extreme cases: in 
one case the disk is marginally stable and spiral waves develop 
(models ml and m2), with the main global effect of heating in the 
radial direction (the Araki parameter (Tz/crr ~ 0.5 at two scale 
lengths from the center at the end of the simulations and it de- 
creases at larger radii, see the middle panel of Fig. [T]). In the other 
case (model m6), the disk is heated and the velocity dispersions 



Table 1. Initial configuration of the A/^-body simulations: name 
of the model, initial Safronov-Toomre parameter at two scale 
lengths from the center, ratio between halo and disk-bulge 
masses, initial radial and vertical velocity dispersions at two 
scale lengths from the center, ratio between hot particles and 
total visible particles. 



Model 


Qt 


Mh/Mbd 


(Tr [km/s] 


<Tz [km/s] 


NhotIN 


ml 


1 




20 


18 


12% 


m2 


1.5 




30 


20 


15% 


mS 


2 


2 


30 


20 


4% 


m4 


4 


4 


30 


22 


2.5% 


m5 


5 


4 


40 


22 


2.5% 


m6 


5 


4 


40 


40 


5% 



increase in both the radial and the vertical direction (the Araki 
parameter cr^/cTr ~ 0.8 at large radii at the end of the simula- 
tion, see the corresponding curve in the middle panel of Fig. [T]). 
The other models considered have intermediate values of Qt and 
(TzIcTr between these two extreme cases. The Safronov-Toomre 
and Araki parameters at the final times are shown in the left and 
middle panels, respectively, of Fig.[T] 

In the right panel of Fig. [T] the rotation curves for all the 
A/^-body models listed in Table 1 are shown at the final time. 
When the halo component is included (in models m3 , m4 , m5 , 
m6), the rotation curve is flatter and Vc ~ 220 km s"^ at two 
scale lengths from the center, that is in the region between 6 and 
10 kpc, depending on the model. The rotation curve of model m3, 
where the ratio between halo and disk-bulge mass is Mh/Mbd = 
2, has intermediate values between that of models without the 
halo (ml , m2) and the others. 

We classify stellar orbits into three dynamical categories 
(Sparke & SellwoodTW; Pfenniger & Friedli [1991]). The first 
two dynamical categories are the bar and the disk orbits with 
the Jacobi integral H = E - flpL^ smaller than the value at the 
Lagrangian points Li 2,H< H(Li 2), where E is the total energy 
and Lz is the z-component of the angular momentum. The sep- 
aration of particles in the bar or disk component can be easily 
done since bar orbits typically have smaller values of and E 
than disk orbits. The third category includes hot orbits for which 
H > H(Li 2). The models considered here differ in the number 
of hot particles after the formation of the bar (the ratio between 
the number of stars in the hot component and the total number 
of the visible stars is listed in the last column of Table [TJ. 

The bar pattern speed flp = ^(t), where 6 is the azimuthal 
angle of the bar major axis (in the inertial frame) calculated by 
diagonalising the moment of inertia tensor of the bar particles, 
ranges from 35 km s"^ kpc~^ for model ml to 40 km s"^ kpc~^ 
for m6, at final times. These values are comparable to those found 
by Fux (119971 - see his Fig. 5), while recent estimates of the 
MW bar pattern speed are of the order of 50 - 60 km s"^ kpc"^ 
(Dehnen 2000, Minchev et al. 2007, 2010). The pattern speed 
is typically slowly decreasing in time with a rate of a few 
km/s/kpc/Gyr (Fux 1997. Bournaud & Combes'2002\, so that the 
values at the beginning of the simulations are 60 km s"^ kpc"^ 
and 80 km s"^ kpc~^ for models ml and m6, respectively. The 
corresponding corotation radius Rc, obtained by the intersection 
of with the circular frequency, ^p(t) = ^(Rc, t), increases in 
time (it typically ranges from Rc = 2 to 5 kpc at final times in 
our models). 

In order to understand the role played by the central bar on 
the distribution of stars in the disk, we follow the evolution of the 
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Fig. 2. Evolution in time of the bar's strength for models ml and m6. 



strength of the bar in time for each model. If Cm is the amplitude 
of the mode m in the density distribution, 



r 



^ Qxp(im6j) 



(3) 



the bar's strength is defined as the mode C2 when the stars j are 
restricted to the bar component. We normalise C2 with respect 
to the number of stars in the bar component, Cq. This quantity 
is shown in Fig. [2] for models ml and m6. Since we do not in- 
clude the gas component in our models, we are not able to fol- 
low their evolution for times longer than a few Gyrs. After this 
typical time-scale, the bar amplitude saturates and the systems 
reach a quasi- steady state. Since the inclusion of the gas compo- 
nent necessarily requires the introduction of other less controlled 
parameters, such as the cooling rate or star formation, we prefer 
in this first study to limit the integration time over a couple of 
Gyrs, which is already enough to study the role of the bar in the 
radial migration process. 

The face-on and edge-on views of the density distribution of 
models ml (upper panels) and m6 (bottom pannels) are shown 
in Fig. [3] at time t ~ 550 Myr, that is just after the maximum 
strength of the bar (cf. Fig. O. In the first case, both bar and 
spiral arms develop since the disk is sufficiently cold, while in 
the second case the disk is hot and only a bar (with a smaller 
corotation radius than in the previous case) develops in the cen- 
tral region. In the external regions of models ml and m6, other 
patterns can be observed, the dominant being the pattern with 
m = 1. This mode is related to asymmetric distributions of mass 
pushed by the system rotation and it can give rise to filaments of 
stars or ring-like structures on long time- scales. 



3. Diffusion equation in axisymmetric systems 

In this paper we model the diff'usion along the radial direction in 
the galactic plane by introducing a distribution function F(R, t) 
which satisfies the phenomenological spatial diff'usion equation 
(from Fourier's law) in cylindrical coordinates: 



dtF = ^dR(RDdRF), 
K 



(4) 



where D is the diff'usion coefficient (in unit of area per time). 
This diff'usion model has already been more succinctly presented 
in Brunetti et al. (120101) . We suppose that for a limited time all 
the stochastic processes ongoing in a spiral galaxy can be de- 
scribed by such an equation characterized locally by a single 
positive diff'usion coefficient Z), to be empirically measured in 
the full dynamical simulations for a range of R and t. The pa- 
rameter D conveys the notion that the diff'usion process is a sur- 
face increase per unit time. For the sake of simplicity we choose 
D to be constant and independent of F, which makes Eq. ^ a 
linear partial diff'erential equation for F. Another choice could 
have been to take C = Dp as a local parameter, where p is the 
mass density. The unit of C is mass per length per time which 
conveys then the notion of a mass flux gradient. However this 
choice would have led to a more difficult empirical determina- 
tion of C since the equation involves then p and its gradient. 
Thus, the model considered here gives a lower limit on the dif- 
fusion happening in a barred galaxy, since other more complex 
diff'usion processes are neglected in the present analysis. 

The general solution of Eq. dU with constant D, which is 
non- singular at 7? = is given by: 



-f 

Jo 



F(R, t) = I A(s) e-^'' Jo(sR) s ds , 



(5) 



where /o(-^) = '^o(--^) is the Bessel function of the first kind. The 
function A can be determined by taking the Hankel transform of 
F{R. 0): 



Jo 



A{s)= F(R,0)Jo(sR)RdR 



(6) 



By inserting Eq. © into Eq. (O and assuming that the parti- 
cles are initially localized at a certain radius Rq at time to = 0, 
F(R, 0) = FoRoS(R - Rq), we obtain: 



F(Rj) = RlFo 

Jo 



se-^'' Jo(sR)Jo(sRo)ds 



(7) 



Thus, the diff'usion of this distribution can be expressed in terms 
of Bessel and elementary functions (Gradsteyn & Rvzhik i20071 
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formula 6.633.2). The time-evolution of an initial set of localized 
particles reads: 



2Dt 



exp - 



ADt 



/o(» 



2Dtl 



(8) 



where /o(x) is the modified Bessel function of the first kind, 
which is finite at the origin, /o(0) = 1. In Fig. [4] two initial dis- 
tributions with D - \ and centered in Rq = 0.5 and 2 (solid 
lines) evolve in time, as described by Eq. ([8]) (dashed and dotted 
lines, respectively). At large radii, the distributions are essen- 
tially Gaussian, while at small radii they are strongly modified 
from the contributions of particles at the center of the cylinder. 
Eq. ([8]) describes the distribution of the radial positions of stars 
in the disk at the initial time ti which diff'use toward position Rq 
at time to, such that - ti = At < To, where To is the diff'u- 
sion time- scale or, equivalently, the distribution of stars which 
initially are in Rq at and then diff'use toward R with At <Td- 
When R goes to zero, Eq. ([8]) reduces to: 



R F 

n0j) = ^cxp[-Rl/(4Dt)] 



2Dt 



(9) 



For large values of the argument the modified Bessel function 
Io(x) (2nx)~^^^ exp(x) and thus F(R, t) reduces to: 



lim F{R, t) 

R^oo 



^[AnDtR 
^^R 



exp 



(^-^o)^ 
4Dt 



(10) 



where N(jd, cr) is the Gaussian distribution with mean value fi = 
(R) = Rq and standard deviation cr = ^2Dt. 




Fig. 4. Two distributions with D - \ and centered in Rq - 0.5 and 2 
(solid lines) evolve in time (dashed and dotted lines, respectively), as 
described by Eq. ([8]). 



In order to obtain a simple model of the distribution of the 
stars in a galactic disk, one can consider to envelop Eq. (fTOb 
by an exponential surface density S(7?) oc exp(-R/Rd) (see, for 
example, Sellwood & Binnev .2002) , where Rd is the disk scale 
length, thus obtaining: 



Pd(Rj) 



i47:DtR 
i47:DtR 



C exp 



(R-Ro? 



exp 



(R 



ADt 
■Ro + (T^/Rd? 



Qxp(-R/Rd) 



ADt 
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= C J-^N(M',cr) (11) 

where C and C are normalization constants and N(jd\ cr) is the 
Gaussian distribution with mean value ju' = (R) = Rq - cr^/Rd 
and standard deviation cr = ^j2Dt. It is important to remember 
that the diffusion model used for obtaining Eq. (fTOb is only valid 
for times smaller than the diffusion time- scale. 

We set Rq = Rq = S kpc and Rd ^ kpc. We consider the 
distribution of stars in 7? ~ 7?© at the initial time. The previous 
expression, Eq. (fTTI) , can be used to estimate the relative fraction 
of stars which remain in a diffusion time- scale within the local 
volume |7? - 7?o| <d. This is given by: 

r-Ro+d 

p(\R-RQ\<d)= N(fi\o-)dR (12) 

jRQ-d 

The integral can be written as: 

p(\R-Re\<d)=-^ r%-^' Jx= i[erffc)-erf(x_)] (13) 

where x+ = (cr/Rd ± d/cr)/ V2. If d ^ cr, we get x+ ~ X- and 
thus the probability of staying in the local volume is nearly zero. 
If J ~ o- « Rd, we have p(\R - Rq\ < J) ~ |[erf(l/ ^^2) - 
erf(-l/V2)] = erf(l/V2) = 0.68. If J ~ cr ~ 7?^, we have 
p(\R - Rq\ < J) ~ |erf(V2) ~ 0.48. Thus, the fraction of 
stars which remain in the local volume in a diffusion time- scale 
strongly depends on the ratio d/cr and on the value of Rd. 

The diffusion coefficient D in Eqs. (O-® can be regarded 
as an instantaneous coefficient which depends on the position at 
which particles are initially localised and on the diffusion time, 
since, as already mentioned, modeling the stellar migration as 
a diffusion process is valid only for time intervals less than the 
diffusion time-scale. At <Td- As we described in Brunetti et al. 
(120101) , the diffusion time- scale can be estimated from the simu- 
lation results and it turns out to be of the same order of the rota- 
tion period, Td ~ Trot = 27r/Q(7?, 0- The diffusion coefficient is 
calculated by applying the nonlinear least- square method which 
minimizes the difference between the numerical results and the 
general solution of the diffusion equation described by Eq. ([8]) 
for times less than To- At each time and radial position in our 
N-body simulations, we estimate the instantaneous diffusion co- 
efficient and related quantities, thus obtaining a description of 
stellar migration along the whole simulation. 

4. Results 

The equations and the numerical methods described in the previ- 
ous section and in Brunetti et al. (12010b allow us to calculate the 
diffusion coefficient D(R, t) which is shown in the contour maps 
of Fig. [5l top row, for the models ml (left panel) and m6 (right 
panel). The others models m2, . . . , m5 have intermediate values 
of the diffusion coefficient between these two extreme cases. The 
bar's corotation radius is shown in Fig. [5] as a dashed white line. 
It can be seen that the diffusion coefficient is not constant in time 
nor in radius. If the disk is not too hot, D has the largest values 
D ~ 0.12 kpc^Myr"^ outside the corotation radius of the bar 
(which increases in time in our simulations), where the density 
is strongly perturbed by a m = 2 pattern created by the bar and 
the transient spiral arms, and in the external regions 7? > 8 kpc, 
where the density is modulated by a m = 1 pattern. The stars 
respond collectively to these modulations and the process of mi- 
gration corresponds to a diffusion in an axisymmetric system. In 



hot disks, the diffusion coefficient is large only in the external re- 
gion 7? > 10 kpc, where the m = 1 mode appears, with values of 
the order of Z) ~ 0.08 kpc^ Myr"^ In this latter case, the disk is 
not sufficiently cold to respond to the m = 2 perturbation created 
by the central bar. 

When the disk is marginally stable with initial Safronov- 
Toomre parameter 2r ~ 1, the diffusion coefficient has the 
largest values just outside the corotation region, where two dif- 
ferent families of orbits are present, as can be inferred by the 
total energy and angular momentum values of the particles in 
this region. We consider three different bins of particles located 
respectively in 7? = (1.5 ± 0.5) kpc, R = (3.0 ± 0.5) kpc and 
R = (8.0+0.5) kpc at t = 2.2 Gyr in the model ml. In the first bin, 
particles are mainly inside the bar. Their energies E span small 
negative values and the z-component of the angular momentum 
is nearly centered in ~ (see the two panels of Fig. [6l blue 
lines, labeled as 'bin 1'). Particles in the second bin centered in 
R = (3.0 + 0.5) kpc at t = 2.2 Gyr belong to two different types 
of orbits: one family can migrate only inside the bar, the other 
can go outside the bar, in the disk. Large values of the diffusion 
coefficient D near the corotation region (see the left panel of 
Fig. [21 top row) are related to this superposition of two families 
of orbits. The corresponding values of E and are shown in the 
panels of Fig. [6l (red lines, labeled as 'bin 2'). The two peaks in 
the energy distribution correspond to the two families of orbits: 
bar particles have large negative energies, while particles which 
can go into the disk have small negative energies. Intermediate 
values correspond to the so called hot particles, which as already 
mentioned before have a Jacobi integral H larger than the value 
of H at the Lagrangian points L12, H > H(Li 2). We will dis- 
cuss the hot particles later on in this section. It is important to 
note that the two peaks are increasingly less evident in hotter 
disks and the number of hot particles decreases as Qt increases 
(see the last column in Table 1). In the third bin, particles belong 
to the disk component: these disk particles have small negative 
energies and large values of (black lines, labeled as 'bin 3'). 

From D we can estimate the radial dispersion cr in a rota- 
tion period (which is of the order of the diffusion time- scale), 
cr = ^|2DJ\^i. In the middle row of Fig. [5] we show the contour 
maps of the radial dispersion for models ml (left panel) and m6 
(right panel). If the disk is sufficiently cold, the radial dispersion 
is high near the corotation region of the bar, where it can recur- 
rently assume values of the order of cr ~ 6 kpc. This implies that 
internal stars can recurrently be forced by the activity of the bar 
to migrate in the external region of the disk. At an intermediate 
radius, such as 7? = 6 kpc, the radial dispersion is of the order 
of cr ~ 3 kpc and it increases in the external less dense regions 
where the pattern m = I dominates. In the external region, the 
diffusion of the stellar component is related to the presence of 
both patterns m = 1 and m = 2 and it can be enhanced in regions 
where the bar's outer Lindblad resonance overlaps with the spi- 
ral arms resonances (Minchev et al. 1201 I t). 

The radial migration driven by the bar seems to be efficient 
in cold disks. According to our results for model ml, the number 
of stars which stay always in a local volume of 100 pc around 
7? = 8 kpc is low, since J - 100pc^cr~5 kpc (see discussion 
after Eq. ([T3])). If the disk is hot, the values of cr are lower than 
the corresponding values for cold disks (at 7? = 6 kpc the radial 
dispersion is of the order of cr ~ 1 kpc for model m6 (see Fig.[5l 
middle row, right panel) and consequently radial migration is 
less effective than before. 

Another interesting quantity is the diffusion velocity in a ro- 
tation period, defined by vd = cr/T^ot = ^f2DJJ\^l. We show the 
contour maps of this quantity in the bottom panels of Fig. [5] for 
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Fig. 5. Top row: Contour maps of the diffusion coefficient D. Middle row: Contour maps of the radial dispersion cr = ^|2DTrot^ Bottom row: 
Contour maps of the diffusion velocity Vd - crlTyot - ^/2D/Trof Left: model ml, right: model m6. 



models ml (left) and m6 (right). The diffusion velocity can reach 
values of vd ~ 40 km/s near the corotation region in the model 
ml, while it is always less than vd ~ 20 km/s in the model m6. 

The fact that the corotation region plays a crucial role in stel- 
lar diffusion can be also seen from Fig. El Here it is shown how 
the radial position of stars at some final time, 7?now at tnow = 
2.2 Gyr (on the vertical axis), depends on the corresponding ra- 
dial distribution 7?past at time tpast = 300 Myr (on the horizontal 
axis). The different colors in the map are related to the ratio of 



the number of stars at the two times, N(R 



past-) I- past 



)/mR, 



now-) ''now 



It can be seen that in the case of model ml (left panel) particles 
near the corotation radius, 7?now ~ = 4 kpc (which is the 
value of the corotation radius at tnow = 2.2 Gyr in the considered 
model) came from inside and outside the corotation region, since 
they belong to two different families of orbits, as discussed be- 
fore (see Figs. [6]). Particles outside the corotation at the final time 
were spread over the disk in the past, with Rc < Rp^st ^10 kpc, 
while particles well inside the corotation radius, Rnow < 2 kpc, 
were confined into the bar also in the past. On the contrary. 



in the hot model m6 (see the right panel of Fig. O particles 
were essentially located at the same positions in the past, with 
^past ~ ^now ± AT?, where AT? ~ 1 kpc and it slightly increases 
with Rnow In this case, no particular activity is observed in the 
corotation region, T?c ~ 2 kpc. 

As a further example, let us consider a bin of stars located in 
f^now = (8.0 + 0.1) kpc at time tnow = 2.2 Gyr. Their evolution 
history is shown in Fig. O for the two models, ml (left panel) 
and m6 (middle panel). When the bar's strength is maximum (at 
t ~ 350 Myr, cf. Fig. [2]), stars localized near the corotation radius 
are forced to move toward larger radii, since there the diffusion 
coefficient is large (cf. Fig.O top row). Most of the stars bounce 
then back and forth in the region Rc < R < 10 kpc, to reach 
the final bin position. The evolution history in heated disks such 
as that of model m6 (see Fig. [8l middle panel) is very different, 
since stars are always localized in the region R = Rnow ± cr, with 
Rnow = 8 kpc and cr ~ 2 kpc, in agreement with the correspond- 
ing value of the radial dispersion which can be inferred from the 
right panel of Fig. [5] (middle row). 
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without hot particles. 



Hot particles, characterized by a Jacobi integral H > H(Li 2), 
have a distribution which is maximum in the corotation region. 
In order to investigate the role of such hot particles in the dif- 
fusion process, we have removed them from the bin of stars lo- 
calized in Know = (8.0 ± 0.1) kpc at time tnow = 2.2 Gyr which 
we have considered just before. In the right panel of Fig. [8l the 
evolution history of the bin where the hot particles have been 
eliminated is shown. It must be compared with the evolution of 
the bin which includes them, shown in the left panel of the same 
figure. It can be seen that the main diff'erence between the two 
panels is in the relative number of stars which are able to reach 
the corotation region. The hot particles migrate radially much 
more than the other particles in the disk. The number of hot par- 
ticles in cold disks (such as ml) is nearly three times larger than 
the number in hot disks (see last column in Table [T]). 



2.2 Gyr are in R^, 



B,[kpc] 

■ (8.0 + 0.1) kpc: (left) ml, (middle) m6, (right) ml 



5. Discussion and conclusions 

From the previous sections it is clear that the amount of radial 
migration in disk galaxies is strongly dependent on the bar and 
spiral strengths. As we analyzed N-body simulations without 
gas, in order to focus our study on the pure eff'ect of the bar, 
we can only trace radial migration for -2-3 Gyr. It is out of the 
scope of the present study to compare our simulations with the 
observations of the Solar neighborhood, which are the result of 
10 Gyr of evolution. Indeed, our models were not intended to 
reproduce the conditions in our Galaxy, thus the precise values 
of, for example, the diff'usion time-scale To and the radial disper- 
sion cr due to radial migration obtained from these models do not 
correspond to those valid for the Milky Way. However, the phys- 
ical processes at play in cold disks, like our model ml, should be 
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similar to those happening in our Galaxy, where the Safronov- 
Toomre-parameter in the Solar vicinity is Qt ~ 2. Thus, we 
expect that the order of magnitude of the relevant quantities ob- 
tained from our numerical results are reasonably valid also for 
the earliest phases of the thin disk of our Galaxy (within ~ 2- 
3 Gyr from the bar formation). 

5.1. Implications for chemical evolution models of the Milky 
Way 

Chemical evolution models of our Galaxy traditionally assume 
that the majority of stars do not migrate over large distances, 
and model the Galaxy by introducing independent radial an- 
nuli which are wide enough (around 1-2 kpc wide) so that 
this approximation would be a valid one (van den Bergh 19621 
Sc hmidt .1963; Pagel 1997. Chiappini et al. 1997, Chiappini et 
al. 12001b . The expectation is that intruders from other galacto- 
centric distances would not represent more than a few percent of 
the stars in the local samples. However, as discussed in Sect. 1, 
there are recent claims that radial migration was more efficient 
than previously assumed. This, in turn, is driven by the large 
scatter in the AMR of the Geneva-Copenhagen sample. 

In the particular model of Chiappini et al. (12001b , each annu- 
lus is 2 kpc wide (i.e. d = 2 kpc). Thus, in this case in the Solar 
vicinity d is nearly half the value of the radial dispersion ob- 
tained in model ml at intermediate radii 6-8 kpc, d ~ 2cr. Thus, 
the relation between the different length- scales which is approx- 
imately valid in the Solar vicinity is J ~ cr ~ (see discussion 
after Eq. ([13])) and the percentage of stars which stayed in a vol- 
ume \R - Rq\ < d = 2 kpc turns out to be of the order of 50% 
in a diffusion time-scale. To ~ T^ot = 27iRq/Vc ~ 223 Myr near 
the Sun. The region from which the rest of stars comes from de- 
pends on the activity of the bar which is not a constant pattern, 
as assumed in the past (Wielen 1977 ). If we consider for exam- 
ple the recurrent bar scenario described in Bournaud & Combes 
(12002b , the bar can be rebuilt several times in a Hubble time in 
galaxies with significant gas accretion. When the bar strength 
is high, such as at the beginning of our simulations and in gen- 
eral when the bar is rebuilt by episodes of dissipative infall of 
gas, stars come mainly from the corotation region (see Fig. [8l 
left panel, t < 500 Myr), which is closer to the center since the 
corotation radius typically becomes smaller after each reforma- 
tion episode (Bournaud & Combes 2002 ). When the bar strength 
saturates to a constant value during a quiescent phase, stars can 
span the region between the corotation radius and -10-11 kpc 
from the galactic center (see Fig.[8l left panel, t > 500 Myr). 

In Fig. [9] we show two examples of stellar diffusion. Stars 
localized in 7? = (8 ± 1) kpc and R = (3 ± I) kpc at the end of 
the simulations are in the black bins, in the left and right panel, 
respectively. The red and blue distributions correspond to the 
radial positions of the stars 2.2 Gyr before, for model ml and 
m6, respectively. As can be seen from the left panel in Fig. [9l 
when the disk is sufficiently cold (model ml, red distribution), 
the radial dispersion is higher than the corresponding values for 
the hot disk (model m6, blue distribution). In 2.2 Gyr (which is 
much larger than the diffusion time-scale), only 25% of the stars 
remain in the black bin, all the others come from outside. When 
the disk is hot, the percentage is nearly twice. In the right panel 
we show the case of stars near the corotation region (remem- 
ber that the corotation radius is Rc = 2 kpc for model ml and 
1 kpc for model m6). Stars come mainly from the corotation ra- 
dius when the disk is cold (red distribution), while more than 
60% stay at the same position if the disk is hot. 



We can conclude that the dynamical effects of stellar mi- 
gration should be included in chemical evolution models of our 
Galaxy insofar as the distance d between each annulus is smaller 
than the radial dispersion cr, which is related to the diffusion co- 
efficient D and to the diffusion time- scale Tohy cr = ^J2DT^. 
The radial dispersion a depends on the degree of marginality of 
the disk. 

5.2. Conclusions 

Disk galaxies are complex systems where collective phenomena 
give rise to the emergence of nonlinear structures, such as central 
bars and spiral arms. We have investigated the role of bars in 
marginally stable disks and overheated disks, and their forcing 
effect on the hot (chaotic) particle component, which is more 
sensitive to external/internal perturbations. 

Modeling the migration of stars in marginally stable disks 
as a diffusion process in the radial direction is a powerful tool 
which allows us to estimate quantitatively two crucial parame- 
ters, the diffusion coefficient and the diffusion time- scale. With 
these quantities we are able to compute two other fundamental 
quantities which are the radial dispersion and the diffusion ve- 
locities at different radii and at different times. It is important 
to note that such a diffusion model makes sense only if there 
is a stochastic microscopic component at the origin of the diffu- 
sion. Ideally, the diffusion time- scale should not be much shorter 
than the microscopic e-folding time due to chaos, which in N- 
body systems is typically of the order of the dynamical time 
(Miller 1964). Over longer time-scales, the diffusion model be- 
comes influenced by the mere global dynamical evolution of the 
disk, so its behavior may depart from a simple linear diffusion 
equation with constant coefficient. Thus, the present diffusion 
calculation is useful for diffusion time- scales in a range around 
the rotational period. At each time and radial position in A/^-body 
simulations, we are able to estimate the instantaneous diffusion 
coefficient and the related quantities (i.e., diffusion time-scale, 
radial dispersion and diffusion velocity), thus obtaining a de- 
scription of the stellar diffusion on the whole simulation. 

We have found that the diffusion time- scale is of the order 
of one rotation period and that the diffusion coefficient D de- 
pends on the evolution history of the disk and on the radial po- 
sition. Larger values D are found in cold disks near the corota- 
tion region, which evolves in time, and in the external region, 
where asymmetric patterns develop. Marginally stable disks, 
with 2r ~ 1 , have two different families of bar orbits with dif- 
ferent values of angular momentum and energy E, which de- 
termine a large diffusion in the corotation region. In hot disks, 
Qt > I, stellar diffusion is much more reduced than in the case 
of marginal disks. 

The calculations of both the diffusion coefficient and the dif- 
fusion time- scale give us a quantitative measure of the migration 
process in the disk. Another advantage of studying the diffusion 
of stars in real space, rather than in velocity space, is that it can 
be more easily related to the evolution of chemical elements, 
which can be modeled as tracers which follow the evolution of 
the stellar component. The diffusion process of stars and tracers 
can be directly implemented in chemical evolution codes, which 
we plan to do in the near future. 

It is interesting to compare our results with those obtained 
in a recent analysis of Shevchenko (12011b , where the Lyapunov 
and the diffusion times are estimated for the Quillen's model 
(Quillen 2003 ) which describes the Hamiltonian motion in the 
Solar neighborhood due to the interactions of bar and spiral arms 
resonances. He found that the Lyapunov time, of the order of 
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10 Galactic years, depends weakly on the model parameters, 
which can radically change the extent of the chaotic domain (as 
we obtain by varying the Safronov-Toomre parameter Qt). The 
diffusion time, which characterizes the transport in the chaotic 
domain of the phase space, is calculated as the inverse of the 
diffusion rate in the energy variable (thus differing from our def- 
inition), with upper bounds of the order of 10 Gyr. It strongly 
depends on the radial position in the Galaxy, in agreement with 
our findings. 

We call attention to the fact that although there are good rea- 
sons (both theoretical and observational) to expect that radial 
process took place during the evolution of our Galaxy, it could 
have been much weaker than what has been proposed so far, as 
implied by the existence of radial abundance gradients, and its 
mild time- variation. Unfortunately, the uncertainties in the ob- 
served abundance gradient evolution are still large and we need 
to wait until better ages and distances will be available, which 
will happen in the near future thanks to GAIA and asteroseis- 
mology. 

Meanwhile, the theoretical work should focus on the rela- 
tive importance of the main drivers of radial migration, and in 
the case of the Milky Way, on the role of the bar. Here we 
have shown that the radial migration process is not only time- 
dependent but also changes with galactocentric distance, in con- 
nection to the bar, plus spiral arms. We plan to analyze simula- 
tions with bar but including also gas accretion, where the radial 
migration process can be traced for several Gyrs. In this way 
we will be able to answer if radial migration could have had re- 
peated peaks along the MW history or if it faded away after the 
first 2-3 Gyr. These models coupled with the chemical informa- 
tion will be essential to interpret the radial mixing effects on the 
local age-metallicity relation, the metallicity distributions of the 
local thick and thin disks, and on the evolution of the abundance 
gradients in both disks. 
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